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Abstract 

We present a Monte Carlo numerical investigation of the Hamiltonian Mean Field 
(HMF) model. We begin by discussing canonical Metropolis Monte Carlo calcula- 
tions, in order to check the caloric curve of the HMF model and study finite size 
effects. In the second part of the paper, we present numerical simulations obtained 
by means of a modified Monte Carlo procedure with the aim to test the stability of 
those states at minimum temperature and zero magnetization (homogeneous Quasi 
Stationary States) which exist in the condensed phase of the model just below the 
critical point. For energy densities smaller than the limiting value U ~ 0.68, we 
find that these states are unstable confirming a recent result on the Vlasov stability 
analysis applied to the HMF model. 
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1 Introduction 



The Hamiltonian Mean Field (HMF) model, introduced in ref. [1,2] is a classi- 
cal many-body model of fully-coupled rotators, that has recently raised a great 
interest for its connections to nonextensive thermodynamics [3,4,5,6,7,8,9] and 
glassy systems [10,11]. It can be seen as a system of N planar classical inertial 

spins Ti = (cos6i, sinOi) which interact through an infinite-range potential or 
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equivalently as N interacting particles moving on the unit circle [2] . Indicating 
with K the kinetic energy and with V the potential energy, the Hamiltonian 
can be written as 

N 2 1 N 

H = K + V = J2Y + w^ [1 - COs{6i - d > )] ' (1) 

i=l i,j=l 



where — n < 9i < tt is the ith angle and pi the corresponding conjugate 
momentum. The degree of clustering of the system can be expressed through 
the usual order parameter M, defined as 

M = ^|E^| • (2) 



i=i 



The model can be solved exactly in the canonical ensemble formalism and ex- 
hibits a second-order phase transition from a low-energy condensed (ferromag- 
netic) phase with magnetization M ^ 0, to a high-energy homogeneous one 
(paramagnetic), with M = 0. The dependence of the temperature T = 2< / ^" > 
on the energy density U = E/N at equilibrium is given by the following caloric 
curve [1,2] 

f/=| + i(l-M 2 ) . (3) 



The critical point is at energy density U c = 0.75 and the corresponding crit- 
ical temperature is T c = 0.5. The dynamics of HMF can be investigated by 
numerical integration of the equations of motion at constant energy. Starting 
the system with out-of-equilibrium initial conditions, for example adopting 
the so-called Ml initial conditions (i.e. considering 9i = for all i - so that 
M(t = 0) = 1 - and velocities uniformly distributed), the results of the mi- 
crocanonical molecular dynamics (MD) simulations show a disagreement with 
the canonical prediction in a region of energy values 0.5 < U < U c [3,5,8]. 
Here, for a transient regime whose length depends on the size N, the system 
remains trapped in anomalous metastable quasi-stationary states (QSS) at a 
temperature lower then the canonical equilibrium one. After such a transient, 
for a finite size N, the system slowly relaxes towards Boltzmann-Gibbs (BG) 
equilibrium, showing aging and power-law correlations [8]. Moreover, in the 
thermodynamic limit, if the infinite time limit is considered after iV — > oo, the 
quasi- stationary states become stable and the QSS regime can be considered 
as a true non-canonical equilibrium phase of the model [3]. In the last years 
many investigations have been performed in order to explain the nature of 
the dynamically-generated anomalies of the QSS regime. The most promis- 
ing scenarios seem to be Tsallis nonextensive statistical mechanics framework 
[3,5,8] and the glassy-like weak ergodicity breaking description [10,11]. Both 
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of these scenarios indicate that, during the QSS regime, the system remains 
trapped in a very complex (fractal) region of the phase space, which hinders 
the complete visit of the whole a-priori accessible phase space. 
In this paper we present a Monte Carlo study of the HMF model. In fact an 
extensive Monte Carlo investigation is missing in the literature and only in 
refs. [12,13] some calculations for very small sizes were discussed. The present 
paper is divided into two sections. In the first one, by means of a standard 
Metropolis algorithm, we reproduce the canonical equilibrium caloric curve of 
the model and study finite size effects close to the critical point where fluctua- 
tions are larger and simulations more delicate. In the second section, we mod- 
ify the standard algorithm in order to perform a sampling over the constant 
energy hypersurface of phase space, looking for those states with minimum 
temperature (kinetic energy), not necessarily Boltzmann-Gibbs equilibrium 
states. Below the critical energy U c = 0.75 density, these spatially homoge- 
neous (M = 0) states coincides with the microcanonical non-equilibrium QSS 
for energy greater than U ~ 0.68. With our modified Monte Carlo algorithm 
we found that below this limiting value, as also showed in a recent paper [15] 
by means of a nonlinear stability test, the homogeneous states are effectively 
unstable. We discuss these results comparing them with the out-of-equilibrium 
caloric curve calculated using molecular dynamics and with the absolute min- 
imum temperature curve. 



2 Metropolis Monte Carlo simulations 

The Metropolis Monte Carlo algorithm is a well-known general method [14] for 
computing the canonical equilibrium statistical expectation values by means 
of a weigthed random sampling of the possible microstates. The algorithm 
usually generates a Markov chain of configurations, for which the probability 
of having a given configuration C n depends only on C n _i and not on the 
previous history of the system. Given a configuration C n , one extracts a new 
trial configuration C with a random algorithm characterized by a simmetric 
transition probability, in order to satisfy the detailed balance condition and to 
minimize the free energy. For example, in the Ising model the configuration C 
could be the configuration C n in which a random chosen spin has been flipped. 
More generally, the configuration C could be obtained from C n by changing 
at random the state of the ith spin of the system considered. Then, if the 
equilibrium distribution is exp[—(3H(C)], where H(C) is the Hamiltonian of 
the system and (3 the inverse of the temperature T (which remains constant), 
one computes H(C) — H{C n ) = AH and uses the extraction of a uniformly 
distributed random number r in the interval [0,1] in order to accept the new 
configuration, according to the following rules 

C n+1 =C n if exp[-(3AH]<r 



3 



.Ill III 1 I I MM 1 I I INI 1 I I I III 

i T=0.2 (a) 



ZD 

>-. 

c 
CD 

T3 

>-, 

cn 
±~ 

03 

tz 
LU 



lc+05 
Se+04 
6c+04 
4c+04 
2e+04 




~i — I — | — I — I — | — I — I — | — I — I — | — I — r 



o. is 



0.22 



0.24 



U mc =0.6798 



0.(0 



0.70 



10 10 

Monte Carlo iterations 



O.DO 

U 



(d) 



0.7> 



=3 
CO 

>, 

o 

tz 

CD 
=3 

tr 

CD 









1 U =0.9 

mc 



















0.D2 



Fig. 1. For the case ./V = 1000 we plot in panels (a), (c), (e) the standard Metropolis 
Monte Carlo evolution of the energy density U and the relative histograms, in panels 
(b), (d), (f) for three different values of temperature, i.e. T = 0.2,0.476,0.8, which 
are , respectively, below, near and above the phase transition. 

C n+1 = C if exp[-(3AH}>r , (4) 



By means of this technique one can generate a sequence of configurations from 
which, after an opportune thermalization transient, it is possible to get the 
canonical equilibrium properties of the system. In fact, if N iter is the number 
of iterations of the algorithm, the expectation value of any observable quantity 
A{C) can be calculated using the identity 



lim ^ £; A(C n ) = [ d»(C)A(C), (5) 

L 71=1 



with 

dfi{C) ~ exp[-/3H(C)], J dfi{C) = 1 . (6) 

At variance with usual statistical models, the HMF model has also a kinetic 
term. Thus in order to use the standard Metropolis algorithm and calculate 
the caloric curve of the HMF model, we fixed the parameter (3 = 1/T, in 
this way the kinetic energy K = NT/2 of the system of N rotators is fixed 
as well and remains constant for all the simulation. We started from initial 
conditions with all the angles equal to zero, i.e. M(t = 0) ~ 1. Then , in the 
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Fig. 2. For the case T = 0.476 which corresponds to the well studied energy density 
value U = 0.69, we plot in panels (a), (b) and (c) the standard Metropolis Monte 
Carlo evolution for different N in order to investigate finite size effects. The average 
value U mc is also reported and shown as dashed line. In panel (d) we plot the values 
of the energy density obtained as a function of the system size. A clear convergence 
towards the theoretical prediction Uth is evident together with a decrease of the 
error bar which is given by the corresponding standard deviation extracted from 
the numerical fluctuations around the average value. 



configuration C n , we changed the angle 9k of a given small quantity (choosen 
properly in order to satisfy the detailed balance condition), so that, in the 
new configuration C, we have 9' k = fe +A0fc. We computed the corresponding 
variation AV in the potential energy and since AK = we have AH = AV. 
Following the rule of eq. (4) we finally accepted or not the new configuration. 
After a termalization transient, we started to compute the average energy 
density U = H/N by means of eq.(5). 

Let us discuss in the following the numerical results obtained with this stan- 
dard Monte Carlo algorithm. In Fig. 1, panels (a), (c) and (e), we plot Metropo- 
lis simulations performed for a system of N=1000 rotators and three dif- 
ferent values of temperatures. We considered in particular the values T = 
0.2, 0.476, 0.8 which are, respectively, below, near and above the phase tran- 
sition crtitical temperature T c = 0.5. After a transient time, the simulations 
converge to a plateau and fluctuate around an average value. We report in 
panels (b), (d) and (f) the relative histograms with the average Monte Carlo 
values U mc obtained. As expected, fluctuations are larger for temperatures 
close to the critical one. 
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Fig. 3. The numerical results of the standard Metropolis simulations performed for 
a system of N=1000 rotators (circles) are plotted in comparison with the analytical 
canonical caloric curve, corresponding to the Boltzmann-Gibbs equilibrium (full 
curve) . 

In order to study finite size effects and fluctuations we discuss the evolution 
for different system sizes at a temperature just below T c . We considered the 
case T = 0.476 which should correspond to the well known energy density 
U = 0.69. In Fig. 2 (a), (b) and (c), we show the Monte Carlo evolution for this 
case for different sizes of the system. While fluctuations are quite consistent 
for N = 100, they diminish as expected with the size of the system. Also the 
value to which the simulation finally converges gets closer to the theoretical 
prediction the bigger the size of the system. This is evident in panel (d) of 
the same figure where the values obtained for different numbers of rotators 
together with the relative standard deviations are plotted in comparison with 
the theoretical expected one. 

Finally, in Fig.3, we plot the average energy densities obtained with the Monte 
Carlo simulations for a wide range of temperature values. The Monte Carlo 
points are compared with the analytical canonical curve (full curve) corre- 
sponding to the Boltzmann-Gibbs equilibrium [1,2]. The agreement is very 
good for this energy size (N=1000). As expected at equilibrium, where ergod- 
icity should be satisfied no anamalies are found. A similar results was obtained 
for microcanonical Monte Carlo simulations [12,13]. On the other hand, as dis- 
cussed in refs. [3,4,5,6,8], the origin of the anomalies is in the dynamics which 
induces frustration and nonergodicity for a transient time - before complete 
equilibration - which diverges with N. In order to try to shed more light on 
these metastable states we modified the standard Monte Carlo algorithm as 
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explained in the next section. 



3 Monte Carlo optimization at constant Energy 

3.1 The algorithm 

We mentioned in the introduction that the so-called QSS, i.e. the out-of- 
equilibrium metastable states that emerge from the microcanonical molecular 
dynamics simulations, become stationary in the thermodynamic limit and lie 
at temperature lower than the canonical equilibrium one. More precisely, up to 
U ~ 0.68, the QSS points lie on the extension, in the condensed phase, of the 
high temperature line of the caloric curve (this line is the geometric place of 
states with M=0, maximum potential energy and minimum kinetic one). This 
means that the homogeneous quasi-stationary non-equilibrium states with zero 
magnetization should be stable only up to those limiting values, but not below. 
In a recent paper on the HMF model [15], the authors apply a nonlinear sta- 
bility criterion (a modification of that one originally proposed in ref. [16]) to a 
selected set of spatially homogeneous solutions of the Vlasov equation, which 
describes the continuum limit of the Hamiltonian Mean Field model. Actually 
these solutions are qualitatively very similar to the zero magnetization QSS 
arising from the microcanonical simulations with Ml initial conditions and 
also the results of the stability test found in [15] were consistent with the nu- 
merical evidence of the disappearance of the homogeneous QSS family below 
U ~ 0.68. 

With the aim to study these anamalous states without using molecular dy- 
namics tecniques, but only performing a sampling over the constant energy 
density hypersurface in phase space, we modified the standard MC Metropo- 
lis algorithm adopted in section 2. Our purpose was not that one of finding 
Boltzmann-Gibbs equilibrium configurations, on the other hand we wanted 
to select those out-of-equilibrium configurations which minimize the temper- 
ature (i.e. the kinetic energy) of the system at constant total energy, in order 
to look for homogeneous QSS's close to the critical point and study their 
eventual stability. The problem is an optimization-like one, i.e. using a Monte 
Carlo procedure, one wants to find the states at minimal temperature with 
the constraint of total energy conservation. 

The new algorithm we have adopted obeys the following rules: 

(1) It starts from initial conditions with fixed magnetization and momenta 
uniformly distributed. Then randomly changes the momentum p k in the 
configuration C n so that pk + Ap k = IT; thus, in the new configuration 
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C, one has 

Pi =Pi + (n ~Pi)Si,k with i = l,....,N; 

(2) it computes the corresponding variation in the kinetic energy 

k' = \ jyi + (n 2 - pi)5 hk ] = K + \tf- \ V \ =► 

(3) then the n + 1 configuration is 

C n+1 = C n if exp[-2AK]<r 

C n+1 = C if exp[-2AK]>r (7) 
where r is the usual random number chosen uniformly in the interval 

[o,i]; 

(4) finally, if the new configuration C is accepted, it calculates the new value 
of the angle 6^ by means of the variation of the potential energy AV = 
—AK, knowing that 

V' = V + AV = ^ Y,[l- 008(91 ~ 0',)] (8) 

and 

e' i = e i + (e-e i )8 itk i = i,....,N . 

After some algebra, one finally obtains the following equation for 

a cosQ + b sinQ = c (9) 
where 

N N N 

a = ^2 cosOi b = sinO-i c = cos(&i — 6 k) + NAK . 

i^-k ij^k i^k 

The solutions of eq. (9) are 

= ±arccos R ± (10) 



with 



R± = ac±bD^ _ 1 < i? ±< 1 D 2 = a 2 + b 2 -c 2 >0. (n) 
a 2 + b 2 
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2 Modified Monte Carlo evolution 

U=0.69 N=1000 




to to 1 to 2 to 3 10 4 10 5 

number o1 iterations 



Fig. 4. We plot four examples of temperature evolution for the energy density 
U = 0.69 and N = 1000 obtained with the modified Monte Carlo procedure at 
constant energy. The curves refer to various initial conditions with uniform distri- 
bution of momenta and magnetization M=l, 0.8, 0.6,0, see text. We report also the 
values of temperature for the canonical equilibrium T eq = 0.476 and for the QSS 
regime in the infinite size limit Tqss = 0.38. 

Thus, from eq.(10), one can calculate the new value of the kth angle, 
taking into account the new constraints (11). Then one can repeat the 
algorithm for the new configuration C and so on, until the system reaches 
a stationary state. At this point, the desired expectation values (in this 
case the temperature of the system at fixed energy density) can be cal- 
culated by means of eq.(5). 

3.2 Numerical results 

With this new Monte Carlo algorithm, we calculated a new out-of equilibrium 
caloric curve for the HMF model, which corresponds to those states with min- 
imal temperature at constant total energy density. 

In our simulations we consider N=1000 rotators and different initial conditions 
with uniform distribution of momenta and different initial magnetization vary- 
ing from M = 1 to M = [8] 1 . Then we follow the system evolution until 
a stationary value of temperature has been reached for each energy density 
value considered. In Fig.4 we show for example the temperature evolution in 
phase space for U = 0.69. The plot shows that, using this new optimization 

1 The inital condition M=0 can be adopted only for U > 0.5 
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Fig. 5. The caloric curve obtained with the Monte Carlo optimization procedure (full 
points) is compared with the Boltzmann-Gibbs canonical prediction at equilibrium 
(full curve) and with the metastable QSS caloric curve (open squares) obtained by 
means of the molecular dynamics microcanonical simulations in the QSS regime. 
We report for comparison also the curve corresponding to the a-priori absolute 
minimal values of kinetic energy (dot-dashed). We can distinguish three different 
regions: (I) In the homogeneous phase all the curves coincide. (II) In the energy 
interval 0.68 < U < U c (U c = 0.75) - marked by the two vertical dashed lines - the 
Monte Carlo curve differs from the Boltzman-Gibbs equilibrium prediction, and lies 
on the extension of the high temperature branch as the microcanonical dynamical 
simulations. (Ill) Below U=0.68, while the molecular dynamics curve rejoins again 
the canonical one, the Monte Carlo curve remains at a lower (but not at its absolute 
minimum) temperature value for the whole condensed phase. 

algorithm, the Boltzmann-Gibbs equibrium temperature T eq = 0.476 is not 
a stable minimum, the curve for M = 1 stays there for a while and then 
escapes to reach a deeper minimum. Instead, as expected for this energy den- 
sity, the most stable solution is always the temperature corresponding to the 
homogeneous QSS Tqss = 0.38, also reported in the figure. 

In fig. 5 we plot the caloric curve obtained with this new Monte Carlo op- 
timization approach (full points). The bars reported indicate the size of the 
fluctuations. In comparison we show also the Boltzmann-Gibbs canonical equi- 
librium prediction (full line) and the metastable caloric curve found by means 
of molecular dynamics microcanonical simulations, performed for N= 100000 
with Ml initial conditions, in the QSS regime (open squares) [8]. It is helpful 
to divide the graph in three regions. 
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In region (I), i.e. in the homogeneous high temperature phase, the three curves 
coincide. This means that in this phase the Monte Carlo optimization method 
is able to recover the expected canonical equilibrium curve. This is in agree- 
ment with the fact that here the total energy is essentially kinetic: the potential 
energy is constant and equal to its maximum value V/N=0.5. 
On the contrary, in region (II) - i.e. between U c and U ~ 0.68 - the Monte 
Carlo curve deviates from the canonical prediction, but coincides with the 
results of the molecular dynamics simulations in the QSS regime. As in the 
previous region, here also the Monte Carlo temperature takes its minimum 
allowed value, consistently with the total energy constraint. This fact is in- 
teresting in order to understand the nature of the QSS regime, which here 
corresponds to a stable attractor of the system when the magnetization (and 
the force acting on the single spin) is zero. We have checked that this result 
does not change with the number of spins considered. 

Below the value U = 0.68 - just lower than the value (U = 0.69) where the 
QSS dynamical anomalies and the disagreement with the canonical prediction 
are more evident - the Monte Carlo curve starts to disagree also with the 
molecular dynamics simulations. In fact, in the region (III), while the QSS 
anomalies tend to disappear and the molecular dynamics curve slowly rejoins 
(below U = 0.5) the canonical Boltzmann-Gibbs equilibrium prediction, the 
Monte Carlo curve stays below the other two curves. However the Monte Carlo 
result does not correspond to the absolute minimum value of temperature. In 
Fig. 5 we report, for comparison, the curve (dot-dashed) corresponding to the 
a-priori absolute minimum value of kinetic energy. The Monte Carlo curve 
differs from such a minimal kinetic energy - which corresponds to M = - 
around the limiting value U ~ 0.68. This is a numerical confirmation of the in- 
stability of the homogeneous quasi-stationary states below this value as found 
in the nonlinear stability test of ref. [15]. 

Finally, it is interesting to observe that below the value U = 0.5, where the 
kinetic energy - and therefore the temperature - could also be null, the Monte 
Carlo curve lies roughly in the middle between zero and the equilibrium tem- 
perature values. This result is not fully understood and could be originated 
by a competition between the two unstable attractors in phase space. Further 
investigations in this respect is required. 



4 Conclusions 

In this paper we have presented a Monte Carlo study of the HMF model. 
In the first part of the paper by means of a standard Metropolis procedure, 
we were able to reproduce the canonical caloric curve of the HMF model 
at equilibrium and study finite size effects close to the critical point where 
dynamical anomalies exist in the out-of-equilibrium regime. In the second part 
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of the paper we studied out-of-equilibrium states by means of a Monte Carlo 
optimization technique. To this end, we have modified the standard Metropolis 
algorithm, in order to obtain a temperature minimization at constant energy 
density and look for those states with minimal temperature, which are not 
Boltzmann-Gibbs equilibrium states. However in this way we could study those 
metastable homogeneous Quasi Stationary States found dinamically below 
the phase transition point. We found that these states are stable for energy 
densities greater than U ~ 0.68. For energy densities smaller than this value, 
the zero magnetization quasi- stationary states are not reached and a different 
caloric curve is obtained. These results confirm what found in [15] by means 
of a nonlinear stability analysis where the authors suggest that this fact is 
due to the loose of stability of the spatially homogeneous solutions of Vlasov 
equation. We hope that this work, supplying a non-dynamical point of view 
and adding new numerical results to the study of the HMF model, will be of 
help in sheding further light on the nature of the several intriguing aspects of 
this model. 

We thank C. Anteneodo, F. Baldovin, V. Latora and C. Tsallis for stimulating 
discussions. We would like to dedicate this paper to Constantino Tsallis for 
his 60th birthday wishing him a still very long and fruitful research activity. 
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